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Abstract 


An upwind-biased, point-implicit relaxation algorithm for obtaining 
the numerical solution to the governing equations for three-dimensional , 
viscous, compressible perfect-gas flows is described . The algorithm is 
derived by using a finite-volume formulation in which the inviscid com- 
ponents of flux across cell walls are described with Roe’s averaging 
and Marten’s entropy fix with second-order corrections based on Yee’s 
symmetric total variation- diminishing scheme. Viscous terms are dis- 
cretized by using central differences. The relaxation strategy is well 
suited for computers employing either vector or parallel architectures. 

It is also well suited to the numerical solution of the governing equations 
on unstructured grids. Because of the point-implicit relaxation strategy, 
the algorithm remains stable at large Courant numbers without the ne- 
cessity of solving large, block tridiagonal systems. Convergence rates 
and grid refinement studies are conducted for Mach 5 flow through an 
inlet with a Mfl compression ramp and Mach 14 flow over a 15° ramp. 

Predictions for pressure distributions, surface heating, and aerodynamic 
coefficients compare well with experimental data for Mach 10 flow over 
a blunt body. 

1. Introduction 

The formulation of an algorithm for obtaining the numerical solution of a system of partial 
differential equations may be divided into two tasks. First, the solution domain is discretized 
by using a collection of points or cells. Relations between quantities at neighboring cells are 
defined such that, in the limit as cell size goes to zero (cell number goes to infinity), the 
numerical relations are consistent with the partial differential equations. Second, a relaxation 
algorithm must be defined which will drive the residual of the approximation scheme to zero 
simultaneously at all the cells in the domain. The two tasks are clearly related, but in the 
parlance of computational fluid dynamics, the first task deals primarily with how the physical 
approximations are defined on the “right-hand side,” and the second task deals with how the 
solution variables evolve on the “left-hand side.” Right- and left-hand sides refer to the position 
of the terms relative to the equal sign in the definition of the complete algorithm. 

When considering the governing partial differential equations for compressible viscous flow, 
the numerical approximations to the component terms can be conveniently defined as a function 
of the physically relevant fluxes and stresses acting on cell walls. This treatment constitutes 
a finite-volume formulation of the governing conservation laws. If properly constructed, it 
guarantees a conservative formulation in that no extraneous sources of mass, momentum, or 
energy will be introduced into the system as a result of numerical imbalances of quantities 
passing from one cell to the next. The numerical approximations to fluxes and stresses across 
a cell wall completely define the finite-volume formulation of the physics on the right-hand 
side. Still, there are an unlimited number of ways to numerically approximate these quantities 
on a cell wall. Viscous, dissipative stresses are generally approximated by Using second-order- 
accurate central differences which numerically reflect the zone of dependence of dissipative 
processes due to the random, thermal velocity of molecules. However, inviscid fluxes are defined 
by more complex wave interactions. (In fact, the physical zones of dependence of viscous and 
inviscid terms are interrelated, but the relationship can be ignored for the purposes of defining 
approximation schemes.) Additional physical insight as to how information is propagated in a 
flow field can be used to better model the physics of the flow. 
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Such insight is provided by a class of approximation schemes referred to as “upwind 
differencing algorithms.” These algorithms can be derived for the inviscid terms by considering 
the Riemann problem defined by the end states at cell centers on either side of a cell wall 
or by performing an eigenvalue analysis of the Jacobian of the flux vector at the cell wall 
with respect to the vector of conserved variables (refs. 1 and 2). The flux at a cell wall can 
be split to account for the contributions of waves coming from the left and the right (refs. 3 
and 4). Such algorithms are most easily defined within the context of first-order-accurate 
approximation schemes which are inherently very dissipative. Simple second-order (or higher) 
corrections to the upwind approximations usually cause severe oscillations in the computed 
solution in the vicinity of high gradient regions in the flow. Often these corrections contribute 
a negative artificial dissipation which overwhelms the natural, positive dissipation of the base 
first-order, upwind scheme and leads to catastrophic instabilities. A physical interpretation 
of these difficulties notes that higher order correction difference stencils will, at times, cross 
discontinuities in the flow or will cross the true zone of dependence of the approximated flux. 
Mathematical analyses of this problem have been used to construct high-order corrections to the 
scheme that will compare fluxes (or flux differences) in the vicinity of cell walls and will choose 
a stencil that ensures a stable, nonoscillatory solution. Correction schemes which approach this 
task in varied ways are known as flux-corrected transport (FCT) (refs. 5 and 6), total variation 
diminishing (TVD) (refs. 7 through 12), and essentially nonoscillatory (ENO) (refs. 13 through 
16). Some of these correction schemes may violate the actual upwind zone of dependence but 
still function properly (ref. 10). 

The selection of a differencing scheme for the right-hand side is generally based on knowledge 
of the scheme’s performance on similar problems. The scheme used in the present study is based 
primarily on the earlier work of Roe (ref. 2), Harten (ref. 7), and Yee (ref. 10). The application 
focus of this work is on hypersonic flows over blunt bodies, including the base and near wake, 
such as would occur on aeroassisted orbital transfer vehicles (AOTV) (ref. 17). Roe’s scheme 
appeared to be well suited for capturing the strong bow shock associated with this application 
where a typical free-stream Mach number is greater than 30. Also, anticipated modifications 
to the basic perfect-gas formulation for the case of real gases and finite-rate chemistry were 
considered to be relatively straightforward. Coakley’s application of a very similar scheme to 
the problem of three-dimensional flow provided additional evidence of the method’s capabilities. 
(See ref. 18.) The evolution of the scheme to the form described herein, with emphasis on the 
blunt body applications, is described in references 19 through 22. 

The second task, as mentioned earlier, deals with the formulation of the left-hand side 
for defining the evolution of the dependent variables. For problems in which a time-accurate 
evolution of the dependent variables is required, the options are restricted with regard to the 
kinds of modifications that can be made in the relaxation algorithm for driving the solution to a 
steady state (assuming one exists). Both spatial and temporal accuracies must be maintained. 
However, if only the steady-state solution is required, one is free to evaluate any element 
of the difference stencil at any iteration (pseudo time) level which facilitates the relaxation 
process. This approach was utilized by Allen and Cheng (ref. 23) in the way they suppressed 
the viscous stability limit in a low Reynolds number flow. Graves (ref. 24) treated variables 
at the cell center of interest and its nearest neighbors at the advanced time step (implicitly ) 
in the “partial implicitization” relaxation technique which greatly enhanced the stability of 
his basic algorithm. In problems utilizing structured grids for multidimensional flow, the 
factored implicit schemes and line Gauss-Siedel relaxation methods permit compromises in 
the way the solution is advanced which significantly reduce the total computational effort while 
retaining solution accuracy (ref. 25). The point-implicit relaxation strategy described herein 
makes similar compromises which make the algorithm ideally suited for problems involving 
unstructured grids and/or executing on massively parallel computer architectures. The essence 
of the strategy is to treat the variables at the cell center of interest implicitly at the advanced 
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iteration level and to use the latest available data from neighbor cells in defining the left- 
hand-side numerics. The success of this approach is made possible by the robust stability 
characteristics of the underlying upwind difference scheme. This strategy was motivated by 
some encouraging results obtained by Lombard et al. (ref. 26) using a sweeping strategy for 
one dimension in conjunction with another upwind algorithm. The basic algorithm requires 
only a single level of storage, and numerical experiments show excellent stability characteristics, 
even when working directly with the steady-state equations (i.e., local Courant number much 
greater than 1). A time- accurate version of the present method is made possible by saving 
at least one extra level of conserved variable data and employing an iterative strategy before 
proceeding to the next time level. 


2. Finite-Volume Fundamentals 


The integral form of the conservation laws applied to a single cell in the computational 
domain is written 




n da — 



u? d£l 


( 2 . 1 ) 


In equation (2.1), the first term describes the time rate of change of conserved quantity q in 
the control volume, the second term describes convective and dissipative flux through the cell 
walls, and the third term accounts for sources or sinks of conserved quantities within the control 
volume associated with thermochemical nonequilibrium. The third term is identically zero for 
the perfect-gas flows considered here. The vectors q and f are defined as follows for viscous 
flow of a perfect gas: 
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The finite- volume approximation to equation (2.1) for a rectangularly ordered, structured 
grid is written 


[^L,* + t ri+i ■ * f * ■ 

+ [ij+i • nj + i<Tj + i — fj ■ a j a j\ IK 

+ [4+1 ■ njfc+l<7fc+l - f k ■ ^*: <7 A:] /J =0 


(2.4) 


where «5q = q n+1 - q n and St = t n+1 - t n . The dependent variable q is defined at cell centers. 
The cell volume ft and cell wall area a are functions of the independent variable r = xf+ yj+ zk 
which is defined at cell corners. A shorthand notation for equation (2.4) that will be used 
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throughout this paper is 


It , + ^ 


l=i,j,k 


^l+l a l+l ~ $ 


Note in equations (2.4) and (2.5) that uppercase variables J, J, K, and L denote computational 
coordinates at cell centers and lowercase variables i, j, k, and l denote cell faces or cell Corners. 
For example, CF i>JK refers to the cell wall corresponding to indices I - J, K, and ajj+l ,j+l,fc 
refers to the cell corner corresponding to indices I + ./ + K — j. A schematic of the 

indexing system is found in figure 1. In the shorthand notation of equation (2.5), the lowercase 
variable l is used as a generic index for i, j , or k. This notation is convenient because most of 
the formulations for quantities at cell faces are independent of the coordinate direction. The 
geometric quantities 0, < r, and n are defined in appendix A. 

The formulations of the inviscid and viscous contributions to the overall conservation laws 
can be considered separately for convenience. Therefore, one can express 

f • n = g + h (2.6) 

where g defines the inviscid terms and h defines the viscous terms. The formulation of these 
terms follows in sections 3 and 4. 

3. Formulation of Inviscid Terms 

The inviscid flux vector at cell face l is defined as 


g / = {o/gL,f + 


(v* • a), 


where 


=[?. 


The first term in braces is a second-order accurate distance-weighted interpolation formula 
for g /. The factors a\ and b\ are geometric weighting functions. These functions, defined as 
follows, account for the relative position of the cell wall with respect to the cell centers: 


Dl + d l- 




(3.2a) 


(3.2b) 


D L = 


°7+i + a l 


(3.2c) 


Geometric weighting improves computed profiles through highly stretched grid regions. An 
earlier version of the algorithm (ref. 21) used volume weighting, in which D L = Vt L , Volume 
weighting tends to magnify the contribution of cells near the axis singularity of spherical 
coordinate systems. - - - 

The second term in braces provides the upwind-biased numerical dissipation. It is a first- 
order dissipation when 0 equals 0 and it is a second-order dissipation when 9 equals 1. The 
variables which define both first- and second-order dissipations are discussed in sections 3.1 

and 3.2. 
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3.1. First-Order Numerical Dissipation 

The vector s ; is defined as 


s, = (Vx ■ n) ; R, 1 {q L - q^) (3.3) 

The matrix R ; _1 in equation (3.3) and the matrices R; and A/ in equation (3.1) are related to 
the Jacobian of the inviscid flux vector g with respect to q in the following manner. 

A - §| = RAR- 1 (3.4) 

The matrices R and R 1 are composed of the right and left eigenvectors, respectively, of 
A, and A is a diagonal matrix containing the eigenvalues of A. These matrices are defined 
in appendix B. The matrix |A| is a diagonal matrix containing the absolute values of the 
eigenvalues of A with constraints on the minimum- allowed magnitude of an eigenvalue described 
in section 6.1. 

The term ■ n J is defined in appendix A. It is related to the ratio of a cell wall 

area to cell volume and may be thought of as the inverse of the projected distance between 
cell centers L and L - 1 on a direction normal to cell face l. The variable y is a generic 
computational coordinate running in the direction of increasing generic index L. Given these 
definitions, the vector s / may be thought of as an approximation to the gradient of the vector 
of characteristic variables (Riemann invariants) in the direction normal to cell face l. Note that 
the term ($x • n) in the definition of s; is canceled by the leading factor of the dissipation 

term in equation (3.1). This leading factor serves to rescale contributions to the second-order 
dissipation and is discussed in more detail in section 3.2. 

The upwind-biased nature of the dissipation term may be understood by considering the 
following difference relation that derives from equation (3.4): 


g L,l ~ EL-1,1 « A;(q£ - q £ _i) 


(v* . it), 


R/Ajs/ 


(3.5) 


The elements of A; or its associated factors R;, A;, and R^ 1 must be determined by 
some suitable average of the dependent variables at the adjacent cell centers L and L — 1. 
Equation (3.5) is only approximately satisfied for general averaging schemes because of the 
nonlinear relation between g and q. However, Roe (ref. 2) has introduced an averaging scheme 
which exactly satisfies equation (3.5) for a perfect gas. This scheme is defined in appendix B. 
The first-order formulation of g/ on an equally spaced grid can now be simplified as follows: 

El = 2 [gL,/ + El—1,1 - R-/|Aj|R ( -1 (qi - qz_i)j (First-order, equal spacing) (3.6) 

The eigenvalues of A, as defined in appendix B, are equal to U, U + a, and U — a where U is 
the normal velocity component through cell face l and a is the sound speed. Consequently, if 
the flow through cell face l is supersonic and positive (moving from cell center L— 1 to cell cen- 
ter L), then all the eigenvalues are positive, |A/| = A*, and equation (3.6) reduces to the 
following relation: 

El — 2 [g L,l + EL-1,1 ~ (g L,l ~ gL-l,/)] = gL-i,l (Supersonic, positive) (3.7a) 


In like manner, if the flow through cell face l is supersonic and negative (moving from cell cen- 
ter L to cell center L - 1) then all the eigenvalues are negative, |A/| = —A/, and equation (3.6) 
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reduces to the following relation: 

g i = \ [g L,l + &L-U + (g L,l ~ g L-l,l)] = SL,l (Supersonic, negative) (3.7b) 

Intermediate situations involving subsonic flow split the contributions to the projected Riemann 
invariants from the right or left according to the sign of the associated eigenvalue. In the most 
general cases, when there are unequal cell sizes, eigenvalues close to zero (to be explained 
in section 6), or second-order modifications, the “upwind” approximation to the inviscid flux 
vector is not quite so exact. However, the present formulation retains the strong stability 
characteristic of upwind schemes over a range of numerical tests that challenge the computation. 

3.2. Second-Order Numerical Dissipation 

A total variation-diminishing (TVD) scheme is an algorithm which is guaranteed to prevent 
the introduction of any new maxima or minima in the evolving distribution of a function or 
functions determined by nonlinear scalar hyperbolic conservation laws or constant coefficient 
hyperbolic systems. The scheme can be second- order accurate and is oscillation free across 
computed discontinuities under these conditions (ref. 10). Although the time-dependent Euler 
equations constitute a nonlinear hyperbolic system, TVD formulations for linear hyperbolic 
systems generally produce oscillation- free solutions when applied to this more complex equation 
set. Conditions for TVD are usually satisfied through the introduction of a limiting procedure 
within the algorithm. For example, several classes of limiters are discussed in references 12 
and 27. The limiters used in the present work were introduced by Yee (refs. 10 and 28). 
They involve symmetric functions of gradients in the neighborhood of the cell face, and 
algorithms based on these limiters are referred to as symmetric total variation-diminishing 
(STVD) schemes. 

STVD schemes involve little extra programming work over simple first-order algorithms 
because most of the quantities required in their implementation are already available, as can 
be seen in the subsequent formulations: 

g m m _ m j n moc j ,si) + min mod (s/_i , s /) - s \ (3.8a) 

sf [n = minmod(s/ +1 ,s;,s/_ 1 ) (3.8b) 

s™ n = min mod (2s m , 2s h 2Si_!, ^ (s/_! + s /+1 )) (3.8c) 

where the min mod function returns the argument of smallest absolute magnitude when all 
the arguments are of the same sign or returns 0 if the arguments are of opposite sign. (If all 
arguments are positive, then a positive result is returned. If all arguments are negative, then a 
negative result is returned.) Also note that equations (3.8) are evaluated element by element 
of vector s. 

The STVD limiter does not yield a strictly upwind bias on the formulation of the flux 
vector because the functional form of equations (3.8) allows for both upstream and downstream 
dependencies. Also, the scheme reduces to first order at cell faces where there is a sign change 
in the arguments of the min mod function. The present formulation differs from that of 
reference 28 in that the differences, Rj -1 (qj: — Ql-x)> have been scaled by (V% • n )/. This 
scaling, which is eventually removed by the leading factor of the numerical dissipation term 
in braces, reduces the effects of grid stretching on the argument returned by the min mod 
function. The present treatment has been found to be well suited to the range of numerical 
test problems considered herein. 
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The dissipation term is defined by the product of the leading factor 


— (which is a 


(Vx • n)j 

first-order difference in the distance between adjacent cell centers projected on the normal to 
the common cell face) and the difference in square brackets [sj - 9s™ m ] (which is zero order if 
9 or sf m equals 0 and is first order otherwise). Thus, the dissipation term is proportional to a 
difference squared when 9 = 1 except for the isolated points in the flow described earlier. 


3.3* Point-Implicit Relaxation of Inviscid Terms 

Equation (3.1) can be approximately linearized with respect to <5q£ in the following manner. 
Define 


8 h ~ { 2 g ^ ! + (°* ~ 2 ) g ^- i + ^ g £-u} 

~ { 2(Vy‘ ■ n), [(** ' fl ), R r‘ (< + ‘ -<d-l) -«*?"'•]} (3.9) 

where superscript n refers to the current value at cell center L, superscript n + 1 refers to 
the new value to be computed at cell center L, and superscript * refers to the latest available 
value at neighbor cell L — 1. The notation gjj refers to the inviscid flux through cell face l 
evaluated with the latest available data from cell center L — 1 and the predicted data at cell 
center L. Elements of the matrices R/, A /, and R. 1 are computed by using Roe’s averaging of 
the current data at cell centers L and L- 1. Elements of the vector (s ; *) min are also computed 
by using current data at cell centers L and L - 1. Substitute g£, + for g 2+ 1 and 

q£ 4- 6<\i for q2 +1 in equation (3.9) to obtain 


- 6; + 2 ( A L,l ~ I A/|) 6q L (3.10a) 

where |A/| = R/|A/|R ( -1 . In like manner, one can show that 

6*+l ,L ~ «+l + 2^ L ' l+l + l A i+lD^^R (3.10b) 

The point-implicit discretization of the inviscid part of equation (2.5) can now be expressed by 

b<\L^L 


St 


~ ~ E Zt+l,L a l+l ~ StL°l\ 


(3.11) 


l=i,j,k 


Use the result of equations (3.10) in equation (3.11) and combine terms to obtain 

| I+ 2^ E [( A i,t+1 + l A /+ll)<R+l - ( A L,l ~ |A/|)<7/] 1 ^q L 
\ l=ij,k J 


~ ~n~ E lSl-hl^l +1 ~ gl^l] 
L l=i,j,k 


(3.12) 


An application of Stokes’ theorem to the summation of and A L[+i in equation (3.12) 
shows that 
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(3.13) 


53 [ A w+i°i+i _ <5q l = Y1 [6%l,i+i<7i+i - L,m\ 

l=i,j,k l=i,j,k 

= 51 ^ f i,inv • Pi+i^i+l - n;<T(] = 0 

l=i,j,k 

Therefore, equation (3.12) simplifies to the following relaxation equation: 

( I + j nv } 6 q L = J] [gf+l^+l ~ Sl^l] 

where 

M L,inv = 2 51 [l A /+ll°’/+l + l A il a /] 

4. Formulation of the Viscous Terms ... 

4.1. Evaluation of Shear Stresses and Conduction 

The viscous stresses on a cell wall with unit normal n in the orthogonal directions n, 1, and 
m are given by 


(3.14) 


(3.15) 



X\ j 

fdU 

dV 

aw a 

. 2 ui dU 

(4.1a) 

Tnn — 

Re * 

{dn 

+ ai + 

dm j 

+ 

' Re 9n 

r nl = 

— 1 
Re 

fdV 

v an 

du\ 

+ ai ) 



(4.1b) 

r nm = 

— | 
Re 1 

fdw 
i, an 

dU\ 
+ amj 



(4.1c) 

velocity « 

components and n. 

, 1 , and m are arc lengths in the n, 1 , 

and m 


directions, respectively. The variables /j, and A are the viscosity coefficients. 

The derivative of some quantity / with respect to arc length in the n direction, for example, 
is expressed as 

at at at at 

(4.2) 


df _ df df df 

- = V/ ■ „ = _„ x + 


Expanding the partial of / with respect to x, y, and 2 in terms of the natural coordinates f , 77 , 
and ( yields 

d£ 

dn 


— + — V77 + —VC 
[dt * dti V 


n 


(4.3a) 


Similar expressions, presented below, can be derived for the derivatives of / with respect to 1 
and m: 


d£ 

d\ 

2i 

dm 


—VC + — V 77 + —VC 

Lac 5 dv v ac 

£^VC + — V77 + —vc 

ac ? dr] v ac 


1 


m 


(4.3b) 

(4.3c) 
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Expanding the derivatives in equations (4.1) according to the formula in equations (4.3) and 
combining like terms yields the following relations for shear stresses: 


_ (2/i + A), f dU - 8U - dU - A 
(«f v{ + ^ v ’' + ac vc ) ■ ” 


Re 


nl “ Re 


Ml 


(4.4a) 

(4.4b) 


i[(** + ** + ?*0MS* + S* + S*)- 

{**♦?*♦?*) -Kf* + ** + **H 

(^ ff+ W f ’ ,+ ? : ^) S+ (f t ’ {+ f^’' + ff^) ■»] 0-4c) 

The component of shear stress acting in the s direction ( s being a dummy variable for x, 
y, or z) acting on a cell wall with unit norm n can be expressed 


T ns ^nn n 5 4* T nl^5 4“ Trunks (4.5) 

Substitution of equations (4.4) into (4.5), collecting terms, and simplifying according to 
the geometric identities provided by equations (A14) and (A15) yield the following relation for 
shear stress in the s direction: 


7"n .«? — 


Hi ((dv-' dv du - \ 

= Rj((s? V{+ at V " + 8C V 9 


4- 


- <9u 


Vtj + 


dn 

9C 


z + d JL d A + Wdii + d_udC\ 

dids dr] ds + dC, ds ) 

■ Vc) n 3 


(4.6) 


where v is a dummy variable for u, v, or w corresponding to s = x, y, or z. Observe that there 
is no functional dependence on the choice of tangential directions f and m in equation (4.6). 

A thin-layer approximation in the y coordinate direction (y = f, r), or () simplifies 
equation (4.6) by neglecting derivatives in the other two coordinate directions. Consequently, 


r ' =£L(^v Y . * + Wdx\ \l( 

Re V<9y d\ ds) Re \ 


dn 

dx 



(4.7) 


where n refers to the direction normal to a surface y = Constant, and the superscript 
prime refers to the thin-layer approximation. The viscous terms on the other two coordinate 
surfaces are also neglected in the thin-layer approximation because their contribution to overall 
momentum and energy balance is small. These approximations are valid as long as the boundary 
layer is relatively thin and the y direction is approximately normal to the high gradient region. 

Equation (4.7) can be further simplified through the use of the geometric relations defined 
by equations (A13) and (A15) and Stokes’ relation (ref. 29), A = — ^ to yield 


The heat flux 


, m (du 1 dU 

ILS Re V<9x + 3 <9y n ' 

through a cell wall is expressed as 
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(4.8) 
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(4.9) 
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(4.10) 


The thin-layer approximation to q n is expressed as 


■ ' _ 

Qn Re Pr 



Let l be the index running in the X direction. Recall that lowercase l refers to cell walls and 
uppercase L refers to cell centers. The solution at cell center L requires evaluation of the 
viscous terms on cell walls l and l + 1. The viscous part of f; ■ n, is now written as 


0 


Tnx 




r n y 
7~nz 


L UT na: + VT ny + W T nz + <7n -I 


(4.11) 


In the case of the thin-layer approximation, h; would be a function of r' and q[ v The viscosity, 
(H, is calculated as the average of n L and m+l- Expressions for the metrics to second-order 
accuracy are presented in appendix A. 


4.2. Point-Implicit Relaxation of Viscous Terms 

Derivatives in the X direction are evaluated to second-order accuracy in space as follows 


(g)’ = «i+i - “2 +1 = “t+i - = 

\OX/l+l,L 



(4.12a) 

(4.12b) 


Terms like ^ for x equal to C are evaluated as follows assuming a rectangular ordering of mesh 
points: 

( jp) = T ( u /, j+i, if “ u /,J-l,if + U *I+1,J+1,K - ( 4 - 13 ^ 

\dvJi+i,J t K 4 v 


Note that the derivatives in the directions along the face (i.e., those derivatives neglected in 
the thin-layer approximation) have no functional dependence on the cell center. Therefore, the 
point-implicit treatment of the full Navier-Stokes equations is identical to that of the thin-layer 
Navier- Stokes equations. 

Now, define h / as a function of differences evaluated with currently available data, for 
example’ (g) ; , and define h* L as a function of differences by using predicted values at cell 

center L, for example (%%)* L - These definitions permit the linearization of the viscous terms 

to be expressed as follows: 


h,* L = h, - B />l 8q L (4- 14a) 

h?+l,L = h i+l+B/+l,L«qL ( 4 - 14b ) 
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where 


B LL = - 


dh 


l,L _ 


dq l 


9Kl 

d^L 


*1+1, L = 


_ dh? +hL 


dqL 


dqL 


(4.15a) 

(4.15b) 


Further details of the definition of matrix B / £, and B i+i t L are presented in appendix C. 

The point-implicit implementation of the viscous terms follows the example set in equa- 
tions (3.11) through (3.15). 


^ - E - »,>] (4-16) 

l=i, j,k 

Note that in the case of the thin-layer Navier-Stokes equations, the summation would only 
include one of the i,j, or k directions, depending on the orientation of the computational 
coordinates with the body. Expand equation (4.16) and combine terms to obtain 

vis j = -— ^2 [hf+l^+l - fycq] (4.17) 

where 

^i.vis = ^ [B/ + 1 1 L<T; + 1 + B/^cq] (4.18) 

l=i,j,k 


5. Relaxation Strategies 


The governing relaxation equation for both the inviscid and viscous components of the 
Navier-Stokes equations is obtained by combining the results of equations (3 14) and (4 17) 
Thus, 



<SQl = r L 


(5.1) 


where 


L M£,inv 4" ^4/, v ; s 

= E (2 l A '+ 1 l + B /+ 1 ,l) a l+ 1 + Ql A il + B /,i) 
l = hjjk ' 




and 

T L = -Q- Kft+1 + h J+l) a l+l - (g / + h /)<q] 

L l=i,j,k 

The relaxed value of q£ +1 can now be determined by 


(5.2a) 


(5.2b) 


-l 

*L (5.3) 

Equation (5.3) involves the inversion of a single 5x5 matrix for three-dimensional flow of a 
perfect gas. In practice, Gauss elimination is used to solve for directly from equation (5.1). 
Note that Mj^inv is obtained as the sum of matrices that have all real, positive eigenvalues and 


q n L +1 =ql + 


ot 

J+ Sh Ml 
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that M/, v i K is obtained as the sum of matrices that have all real, nonnegative eigenvalues. Such 
summations will not generally result in a matrix with real, positive eigenvalues. However, the 
Gauss elimination is routinely implemented without pivoting to facilitate code vectorization in 
the numerical tests, and no evidence of ill-posed behavior has been observed. 

The strategy used to drive the right-hand side of equation (5.1) to zero should take advantage 
of the host computer architecture and the physics of the problem. On computers with serial 
architecture, the relaxation of equation (5.1) (with one or more local iterations) proceeds from 
cell to cell in an ordered fashion. Thus, an updated value of q n + nlocal is obtained at cell 
center (1,1,1) followed by (2,1,1), (3,1,1), . . . (1,2,1), . . . until the dependent variables at every 
cell in the domain have been locally iterated nlocal times using the most recent available data 
from neighbor cells. Vectorizable code is far more efficient than scalar code on machines with 
vector architecture. In this case, all the cells in a plane of the domain are iterated nlocal times 
before proceeding to the next plane. • 

Numerical tests indicate that sweeps which run from a wall across the boundary layer to 
the opposite boundary and then back again are the most efficient for the blunt body problem. 
Effects of a perturbation at a wall are felt at the opposite boundary after one sweep. Effects 
of a perturbation at one cell in a plane parallel to the wall require n iterations to be felt by a 
cell whose index differs from the source cell by n. 

The ordering of the sweeps and the use of local iterations may be used to speed convergence, 
but in numerical tests performed to date, they do not affect the final, converged steady-state 
solution. Thus, it should be possible to solve a large number of cells using a massively parallel 
processing computer in which each cell (or small group of cells) is relaxed semi-independently 
of its neighbor cells (cell groups) using its own processor. The expression semi-independently 
means that a cell (cell group) will need updated information from its neighbor cells (cell 
groups), but neither the order that it receives this information nor the lag time it takes for 
this information to arrive is critically important. (There may be some upper limit on allowable 
communication delays on actual parallel systems before convergence is inhibited. This issue 
is beyond the scope of the present paper.) As long as each processor has immediate access 
to some level of information from its neighbors (which could be stored locally), the execution 
stream could proceed uninterrupted in a parallel, asynchronous mode (ref. 30). 



6. Numerical Issues 


6.1. Zero Eigenvalues and Entropy Violations 

Harten (ref. 7) has shown that the numerical dissipation for waves associated with eigen- 
val ues equal to zero is also zero when the inviscid flux is defined as in equation (3.9). Thus, 
incorrect nonphysical weak solutions of the governing conservation laws which violate the sec- 
ond law of thermodynamics (i.e., expansion shocks) are permitted. In the present algorithm 
it has been found that if no provision has been made for zero eigenvalues, the scheme will 
either fail to converge or will converge to a nonunique solution (usually on a coarse grid) that 
is relaxation path dependent. A simple fix for this problem, suggested by Harten (ref. 7), is to 
restrict the minimum value of the eigenvalue magnitude. Therefore, 



( I ^ieqn I ^ 2e) 

(i-^ieqnl < ' 2c) 


(6.1) 


12 


IiILiIm 



and 


|A| 


Ai 0 0 0 0 

0 A 2 0 0 0 

0 0 A 3 0 0 

0 0 0 A 4 0 

0 0 0 0 a 5 


( 6 . 2 ) 


where A j e(?n is an eigenvalue of A and is defined in appendix B. 

The magnitude of e, which may be viewed as a nondimensional velocity, is problem 
dependent. It is also a function of the method used to nondimensionalize the problem. 
Numerical tests on the blunt body problems, in which there are large zones of subsonic flow 
and stagnation points, indicate that e should be restricted by 


0.1 < e < 0.3 (Approximate range for blunt bodies) (6.3a) 


Numerical tests on flows that are predominantly supersonic indicate that e should be restricted 
by 

0.005 < e < 0.05 (Approximate range for slender bodies) (6.3b) 

The lower limits on e arise from convergence or uniqueness considerations. The upper limits 
arise from concerns regarding excessive numerical dissipation and its effects on predicted surface 
heating, adiabatic wall temperatures, and shock thicknesses. Large values of e have also been 
observed to slow convergence. In some cases for the blunt body problem, the value of e was 
gradually reduced to 0 for cell walls parallel to the body. This treatment improved heat- transfer 
predictions and will be discussed further in section 7. 

Recently, Yee (ref. 12) has suggested a functional dependence of e on the local values of 
sound speed and velocity. This relation has been adapted for use in the present work as follows: 

e^coiai + m + m + m) ( 6 . 4 ) 

where e Q is a constant which generally varies from 0.01 to 0.4. The magnitude of this quantity 
follows similar guidelines as defined in equation (6.3); however, the resultant magnitude of e, 
which is now tied to local velocities, has been found to be generally less dissipative than the 
use of equation (6.3). 


6.2. Relaxation Factors and Frozen Matrices 

Equation (3.14) has been found to be unstable when <9 = 1 (second-order dissipation) 
and the Courant number is greater than 1. (The Courant number, which determines the 
magnitude of 6t , is defined in appendix D.) Because time accuracy is not an issue in the 
present formulation, one is free to adjust M L in any way which improves stability, accelerates 
convergence, or decreases computational time per relaxation step. The instability can be 
removed by multiplying M^ i nv by a relaxation factor $j nv as follows. 

■^L,inv = ^inv^L,inv (6-5) 

where the recommended bounds on the relaxation factor as determined by numerical tests are 
given by 

1-5 < $>inv < 2.0 (6.6) 

Examination of equation (3.14) shows that the application of the relaxation factor reduces the 
magnitude of <5q^ for a given value of the right-hand- side residual. Thus, 4>i nv may be viewed 
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as an underrelaxation factor. It also serves to better approximate the contribution of the flux 
limiter to r^. 

No instabilities have been encountered that are associated with the present treatment of 
the viscous terms. However, an overrelaxation factor associated with the viscous contributions 
to ti can be specified as follows: 


^L,vis — ®vis^®L,vis (6-7) 

where the recommended bounds on the relaxation factor as determined by numerical tests are 
given by 

0.5 < 4> vis (6.8) 

Convergence rates as a function of 4>i nv and $ v j s are presented in section 7, 

Equation (5.2) for is now replaced by 


Mi = M' iiinv + M' i)Vis 
and equation (5.3) for q)' + 1 is replaced by 

q2 +1 =q2-C L r L 


(6-9) 


(6.10) 


where 


C L = 


1+ <t< 


-1 


( 6 . 11 ) 


Note that Cl does not need to be recomputed every relaxation step. In fact, an advantage of 
this formulation is that it requires no more memory or work to save (freeze) C i at every point 
than is required to save M^. This is in contrast to the various classes of implicit methods 
which require the solution of large, block tridiagonal systems of equations in which each of the 
component blocks (which are as large as M^) may be saved, but it is numerically intractable 
to compute and save the inverse matrix. Average computational time per relaxation sweep is 
reduced by a factor of 2 to 3 when is kept frozen for 10 to 40 sweeps. This reduction occurs 
because the work required to recompute inv and vis as well as the work required to 
perform the inversion specified by equation (6.11) is eliminated from approximately 90 percent 
of the relaxation sweeps. The algorithm demands no more work per computational cell than a 
purely explicit formulation except for the effort required to multiply a 5 x 1 vector of residuals 
by a 5 x 5 matrix for three-dimensional flow of a perfect gas. The robust stability for any 
Courant number is not sacrificed with this treatment of C^. 


6.3. Positive-Definite Quantities and Initialization Procedures 

Flow-field initialization is usually achieved by imposing uniform flow conditions at every 
cell. (The only time this initialization procedure failed was in the near wake at cells adjacent 
to the vehicle base in a hypersonic flow. In this case, a linear variation in velocity was 
introduced which varied from 0 at the wall to free-stream values in the middle of the domain.) 
A stability enhancing procedure which has been found useful in the early, highly nonlinear 
solution adjustment period is to force positive-definite quantities (such as density or energy) 
to remain positive definite and to quell any explosive growth caused by physically unrealistic 
initial conditions. The adjustment to the algorithm can be expressed as follows. 
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n+l 
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n+ 1 


= Wl p n 

(6.12a) 

= u>2e n 

(6.12b) 
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where 


i j_ 5qi 

wi = 1 + — 

pTl 

, <5e 

w 2 = 1 + 


<5e 


- [(u 2 + v 2 + w 2 — E) 6 qi — u 8 q 2 — vtfq 3 - w<5q4 + 6 q^ 


us — 


(u> < 0.5) 
(0.5 < w < 2.0) 
(w > 2.0) 


(6.13a) 

(6.13b) 

(6.14) 

(6.15) 


7. Results and Discussion 

The algorithm has been documented in its various stages of development in references 19 
through 22. It is referred to as the Langley Aerothermodynamic Upwind Relaxation Algorithm 
(LAURA). Boundary conditions are discussed in appendix E. Modifications required for an 
unstructured grid are discussed in appendix F. Opportunities for exploiting asynchronous 
relaxation strategies on parallel computers are discussed in appendix G. The material on 
asynchronous iteration and some of the blunt body results first appeared in reference 21. Some 
material on the two-dimensional flow problems appeared first in references 31 and 32. 


7.1. Convergence and Grid Refinement Studies 


Two-dimensional flow problems are used to illustrate convergence properties and solution 
dependence on numerical parameters. The first test case is one studied recently by Rudy 
et al. (ref. 31) involving laminar, supersonic flow through an inlet. The upper wall is straight 
and parallel to the incoming flow. The lower wall is offset 2 cm from the upper wall. Two 
centimeters into the inlet there is a 10° compression followed by a 10° expansion at 4 cm. The 
total length of the domain L is 10 cm. An illustration of the geometry and flow physics in the 
problem is presented in figure 2. The inflow Mach number is 5.0, and the Reynolds number 
based on the length of the domain is 1.14 x 10 7 . Adiabatic, no-slip boundary conditions are 
applied at the upper and lower walls. Extrapolation is used at the outflow boundary. 

Grid refinement studies were implemented for discretizations involving 50 x 50, 100 x 100, 
and 200 x 200 cells. Equation (3.8b) was used to define the limiter, and equation (6.4) was 
used to define e with e Q = 0.1. Cell Reynolds numbers for this test case were evaluated at the 
lower wall in front of the inlet. Values of N c for the three cases are 120, 60, and 30 for the 
three grids, respectively, where 


N c = 


pa Az 
i 1 


(7.1) 


Here, a is the sound speed and Az is the cell size in the direction normal to the wall. This 
definition of cell Reynolds number is used to evaluate the adequacy of grid resolution for the 
purpose of resolving thermal gradients. It depends on a thermal velocity (the local sound speed) 
and not on the local value of a velocity gradient, which can pass through zero at stagnation or 
separation points. Experience with this parameter indicates that it should not exceed a value 
of approximately 2.0 for adequate resolution of the thermal layer. 

The distributions of pressure (fig. 3) and skin friction (fig. 4) show continually sharper 
profiles, particularly in the vicinity of the reflected shock on the upper wall, with increasing 
resolution. It is clear that the separation region in front of the reflected shock (Cf < 0) is not 
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adequately resolved with the 50 x 50 and 100 x 100 cell discretizations. The truncation errors 
cause significant differences in the extent and magnitude of the separation as compared with the 
200 x 200 discretization. The values of N c exceed the recommended limit for adequate resolution 
of dissipative phenomena at the surface. On this basis, the 200 x 200 cell discretization is still 
somewhat coarse. The test cases were repeated with a geometry a factor of 10 smaller than the 
original test. This reduces the free-stream Reynolds number by a factor of 10 and effectively 
reduces the cell Reynolds numbers used above by a factor of 10. (In this case, the finest grid 
is marginally acceptable with a value of N c equal to 3.) These results are presented in figure 5 
for pressure and figure 6 for skin friction. Here again, the finer grids give a sharper pressure 
rise on the upper wall and an earlier separation in front of the reflected shock. The effects 
of truncation error are much less severe (though still observable) for these order 1 values of 
N c (12, 6, and 3, respectively, for the 50 x 50, 100 x 100, and 200 x 200 grids). The small 
differences between the two finest grid solutions can be attributed to the relative coarseness 
of grid transverse to the oblique, reflected shock. Shocks which are oblique to a cell face are 
generally found to be smeared over more cells than in the case of shocks aligned with cell faces. 

A second test case was studied to further address the issue of truncation error. The test case 
involves laminar, hypersonic flow ( Moo = 14.1, Reoo = 2.36 x 10 5 /m) over a 15° ramp. The 
free-stream temperature is 88.9 K (160°R) and the wall temperature is 297.22 K (535°R). This 
case and related tests involving larger ramp angles are extensively discussed in reference 32. 
Four different grids are used in this test: 50 x 50, 100 x 100, 100 x 200, and 200 x 100. The first 
index refers to the number of cells along the wall and the second index refers to the number 
of cells normal to the wall. Values of N c for these four cases, evaluated on the wall in front 
of the ramp at x ss 1.35, are approximately 1.2, 0.6, 0.3, and 0.6. Values of N c near the peak 
skin friction point on the ramp (x « 2.40) are approximately 14.0, 7.0, 3.5, and 7.0. The 
surface pressure and skin friction coefficients as computed on all four grids are presented in 
figures 7 and 8. The boundary layer is well resolved according to the values of N c ahead of the 
ramp and the solutions for both pressure and skin friction coefficients on all four grids are in 
excellent agreement. The solutions on the ramp near the peak pressure and skin friction points 
are in good agreement, though observable differences in magnitude on the order of less than 
5 percent on the three finest grids are observable. The boundary layer thickens but does not 
separate ahead of the ramp. The pressure and skin friction are in excellent agreement among 
all four grids in this region. It is clear that boundary-layer resolution is not the only issue in 
the prediction of surface quantities. The 100 x 100 and 200 x 100 cell cases have identical grids 
in the direction normal to the wall, but differences are observed due to the improved resolution 
across the captured shock in the second case. (The computational speed for the 100 x 100 test 
case is 32 x 10 -6 sec/iteration/grid point on the Cray Y-MP computer using a single processor 
with Cray FORTRAN CFT77 3. 0.2. 2. Approximately 2500 relaxation steps are required for 
convergence.) 

The enhanced numerical dissipation provided by increasing the lower limit on eigenvalue 
magnitude e is shown in the pressure distributions of figure 9 and the skin friction distributions 
in figure 10 in the inlet. These tests were conducted on the 50 x 50 grid for the free-stream 
Reynolds Number equal to 1.14 x 10 7 . The pressure rise across the reflected shock steepens 
with decreasing e 0 . In this coarse grid solution, there is no separation in front of the reflected 
shock; however, the smaller values of e 0 clearly better resolve the sharp decrease in skin friction 
in this area. Computational stability at Courant numbers greater than 5 was compromised 
with values of t Q less than 0.1. Convergence rates deteriorate and solution profiles become 
overly smeared when the eigenvalue limiter is set too large. Nonphysical solutions or a failure 
to converge can result from eigenvalue limiters set too small. 

The effects of the STVD limiters considered herein on the computed solutions are presented 
in figures 11 and 12. The limiters given by equations (3.8a) and (3.8b) differ only where there is 
a local maximum or minimum in the value of s/. The solutions given by these two limiters are 
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nearly identical. The solution given by the third limiter (eq. (3.8c)) gives sharper pressure jumps 
than the first two limiters and lower skin friction predictions. The predictions using this limiter 
are, generally, in better agreement with the predictions given by three other relaxation schemes 
on the coarsest grid, as documented in reference 31 and shown here in figure 13. (The results 
of reference 31 were obtained with a constant eigenvalue limiter equal to 0.005.) Convergence 
characteristics of the scheme with this limiter were not as attractive as those obtained with the 
other two limiters. Usually, the solution would converge to an L 2 norm and then stall at that 
level. Global checks of mass conservation and momentum show gentle oscillations in the fourth 
(or greater) decimal place around a constant level in these quantities. The magnitude of these 
oscillations does not appear to be abated significantly by lowering the Courant number from 
order 10 to order 10 -1 . It would appear that equation (3.8c) is the least dissipative of the 
three limiters based on the sharper profiles in figures 11 and 12 and the oscillatory convergence 
characteristics. 

Convergence histories for two of the test cases are presented in figures 14 and 15 using the 

limiter given by equation (3.8b). The measure of convergence is given by the L 2 norm defined 
by 

1 N 

L2 = 7^Yj T L ■ r L (7.2) 

C N L—\ 

where C N is the Courant number, N is the total number of cells, and r L is the right-hand side 
of equation (5.1). Scaling by C jy reduces the Courant number dependency in the calculation. 
Note that L 2 is not scaled by the total number of cells. Figure 14 shows the rapid convergence 
rate for the coarse 50 x 50 grid solution with a Courant number of 5 x 10® and e u = 0.2. 
Figure 15 shows the convergence history for the fine 200 x 200 grid solution with a Courant 
number of 5.0 and e 0 = 0.1. The minimum value of Cj on the upper wall and the maximum 
value of p / Poo on the upper wall are also plotted as a function of iteration number in figure 15. 
The minimum value of skin friction is found in the separation region in front of the reflected 
shock (fig. 6) and the maximum pressure ratio is found behind the reflected shock (fig. 5). The 
pressure is converged after approximately 600 iterations when the L 2 norm is less than 3 x 10 -3 
and the skin friction is converged after approximately 1200 iterations when the L 2 norm is less 
than 2 x 10 4 in this case. 

An optimum value of $f nv = 1.5 was found from numerical tests on the 100 x 100 grid 
( N c = 6) as shown in figure 16. All these tests were computed with e c = 0.1 and equation (3.8b) 
for the flux limiter. The first curve shows the L 2 norm after 700 iterations from a uniform initial 
flow with 4* vis = 4>j nv and Cjv = 5.0. The solution diverged for values of 4> jnv < 1.4. The second 
curve shows the L 2 norm as a function of 4> inv with 4» vis = 6.0 and C N = 5.0. Tests on other 
problems substantiate the choice of 4> inv = 1.5 as optimum for reducing the L 2 norm. (In some 
tests with particularly severe initialization errors, larger values of 4> iriv would enhance stability 
and convergence during the initial adjustment phase.) The effects of 4> vis and C N on the L 2 
norm after 700 iterations are shown in figure 17. The relaxation factor <h vis has a relatively 
weak influence on convergence over the range tested, with an optimum value occurring at 
approximately 7.0. The large value for $ vis is required as a consequence of stability problems, 
particularly at large Courant numbers, associated with the abrupt transition from free-stream 
conditions at the inflow to no-slip conditions at the wall. The choice of Courant number has 
a stronger impact on convergence, with an optimum value of 5.0 in this test. Optimum values 
for these parameters are expected to be problem dependent, but experience to date has not 
shown significant variation from the values obtained in this test for Courant number or 4> inv . 
Optimum values for $ vis in the blunt body problem discussed in section 7.2 (which has 'no 
boundary condition singularity) vary between 0.5 and 1.0. 
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7.2. Blunt Body Applications on Aeroassist Flight Experiment (AFE) Configuration 

The present algorithm was developed specifically for application to hypersonic flows over 
blunt bodies. The excellent capabilities of total variation-diminishing schemes with regard to 
shock capturing, resolution of severe expansions, and robust stability characteristics motivated 
the development of the present code for these applications. 

The Aeroassist Flight Experiment (AFE) configuration is a 60° elliptic cone raked at 73° 
on the base to form a circular backplane. The nose is defined by an ellipsoid of eccentricity 
equal to 2.0. The shoulder region is defined by circular arcs on planes passing through the 
axis of the elliptic nose. More details of the forebody definition are given in reference 33. An 
illustration of the vehicle is presented in figure 18. The body size is defined relative to the base 
plane diameter. The base plane diameter is equal to 3.9116 m (154 in.) for the flight vehicle. 
The base plane diameter for the wind-tunnel test model is 9.3218 cm (3.67 in.). The inviscid 
relaxation factor $ inv was set equal to 1.5. The viscous relaxation factor $ vis was decreased 
from 1.0 to 0.5, the Courant number was decreased from 5.0 to 0.5, and the eigenvalue limi- 
ter e 0 was decreased from 0.4 to 0.1 over the period of convergence. (There was no comprehen- 
sive retesting of numerical parameters to achieve optimum efficiency.) Free-stream conditions 
are defined by Voo = 1429 m/s, Poo = 60.136 N/m 2 , Too = 52.22 K, and M oo = 9.86. The unit 
Reynolds number is 17040/cm (43 282/in.). Wall temperature is 300 K. 

Predictions from the LAURA code and the HALIS code (ref. 34) are in excellent agreement 
with experimental data for the pitching-moment coefficient in figure 19. The contributions of 
viscous forces to the pitching moment are negligible at this test condition. 

Calculations and experimental data (ref. 35) for pressure coefficient and heat-transfer 
coefficient are presented in figure 20. The data were run at angles of attack of 0°, 5°, and 
-5°. The heat-transfer coefficient for the cone flank compares well with the experimental data; 
however, as angle of attack is decreased from 5° to -5°, the comparison in the stagnation 
region gets progressively worse. At an angle of attack of -5°, where the stagnation point 
crosses the axis singularity, there is an abrupt overshoot/undershoot in the data crossing the 
axis singularity. A grid refinement involving double the number of points across the shock layer 
eliminated most of the irregularity in this region and improved the agreement with experimental 
data in figure 20(c). A grid restructuring, in which the axis singularity was removed by 
redistributing cell faces across the plane of symmetry, also smoothed the prediction curves 
for heating in this region. The strong influence of grid structure on computed heating indicates 
that the variation of truncation error around the axis, as well as the magnitude of the truncation 
error, is responsible for the irregularities in the stagnation region heating. Computed pressure 
distributions are insensitive to these effects. The grid refinement shows little change in the 
heating on the cone flank, which indicates that truncation error is not a significant factor 
elsewhere in the flow field. 

Heat-transfer distributions on the sting in the near wake (ref. 36) compare well with 
calculated data in figure 21. A grid refinement in which a solution was generated with double 
the number of cells in all four domains in the direction normal to the forebody and sting shows 
only a slight change in the computed solution for heat transfer in figure 21(b). The calculation 
assumes laminar flow. There is no direct evidence in the experimental data to indicate if the 
shear layer flow or the sting boundary-layer flow is laminar, transitional, or turbulent. The 
experimental data show no evidence of transition to a turbulent boundary layer on the sting 
(over the length that data were taken) as judged by the heat-transfer distribution. The value 
of the eigenvalue limiter e 0 was set to 0.010. The near wake flow field does not require the 
larger values of eigenvalue limiters needed for stable computation of the stagnation region of 
blunt body flows. 

The four-domain grid used for this case is presented in figure 22 on the symmetry plane 
and on the surface. Grid 1 is the forebody grid which consists of 39 x 24 x 64 cells. Grid 2 
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extends out from the outflow boundary of grid 1 and consists of 31 x 24 x 32 cells which 
approximately span the region extending from above the free shear layer across the captured 
bow shock. Grid 3 fills in the region between the cylindrical upper surface of the afterbody 
with 31 x 24 x 19 cells and spans much of the free shear layer. Grid 4 resolves the boundary- 
layer flow over the sting and the region behind the afterbody with 26 x 24 x 21 cells. The low 
densities in the near wake region simplify the task of constructing grids which satisfy the cell 
Reynolds number requirement at the body surface discussed after equation (7.1). The cell walls 
on adjoining domain boundaries are aligned to simplify flux calculations. The outflow boundary 
condition of the forebody grid system (grid 1) employed second-order accurate extrapolation. 
The inflow boundary conditon for grid 2 picked up from the converged solution in grid 1. This 
boundary condition precludes feedback of information from the wake to the forebody and was 
implemented out of convenience rather than necessity. It is considered reasonable because the 
information from grid 2 can only influence the solution in grid 1 through the subsonic portion 
of the boundary layer which is very thin in these test cases. Grid 2 was constructed so that 
exactly two cell walls of the common boundary in grid 1 formed a single cell wall in grid 2. 

Contours of pressure and Mach number in the symmetry plane at angles of attack equal to 
0° and —5° are shown in figures 23 and 24. Contour lines pass smoothly across grid domain 
boundaries. Gaps in the contour lines are an artifact of the data base supplied to the plotting 
routine. Information is stored at cell centers. The gaps in the contour lines equal the distance 
between cell centers of adjacent domains. The captured shock in the near wake over the sting 
is clearly visible in figure 23. The free shear layer extending from the aft corner of the shoulder 
is evident in the Mach contours in figure 24. 

7.3. Some Observations of the Relaxation Process 

The variation of Courant number, eigenvalue limiter, and viscous relaxation factor noted 
in section 7.2 was implemented to improve convergence. The solution for the AFE forebody is 
initialized assuming that the body materializes in the undisturbed flow at time t = 0. A coarse 
grid, first-order accurate solution is developed on a grid which adapts to the position of the 
captured shock (ref. 21). When the coarse grid solution is judged sufficiently developed, the 
grid is refined and the second-order flux limiters are implemented. A Courant number equal 
to 5.0, a viscous relaxation factor <F v j s > 1.0, and an eigenvalue limiter e a > 0.2 are typically 
applied to speed convergence and maintain sufficient dissipation. Grid adaption is turned off 
once the L 2 norm is observed to be level. Further relaxation can be used to drop the L 2 norm 
about another order of magnitude, but then convergence stalls. At this point, the computed 
surface pressure distribution is generally quite satisfactory, but the computed surface heating 
is likely to differ from the fully converged values by 25 percent or more, depending on the 
amount of grid stretching across the boundary layer and the magnitude of e 0 . The L 2 norm 
can be reduced several more orders of magnitude by reducing the Courant number to 0.5 and the 
viscous relaxation factor to 0.5. The eigenvalue limiter is reduced as much as possible, generally 
to a value of 0.1, in order to minimize the numerical dissipation on the computed heating. The 
eigenvalue limiter can be turned off in the direction crossing the boundary layer in order to 
remove any influence of this parameter on the computed heating. The computation remains 
stable with this adjustment when the solution is close to convergence and it is recommended 
when there is severe grid stretching (for example, when ? BlAK > j 5 w h e re K is the index 
crossing the boundary layer). This adjustment was not used in the present set of calculations. 

The convergence rate in the three wake region domains is better than the rate obtained in 
the forebody region. Typical convergence histories are presented in figures 25 through 30 for 
the forebody and wake regions to illustrate this point. Figure 25 shows the initialization of 
the case for a = 5° on the regular grid starting from a converged case for a = 0° considering 
only the forebody. The Courant number used in this case was equal to 5.0, and e 0 was equal 
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to 0.25. The grid was allowed to adapt to the moving shock during the first 500 relaxation 
steps and then was held fixed for the remainder of the test. Convergence has obviously stalled 
at an L 2 norm of approximately 0.02. The curve in figure 26 picks up from the stalled results 
of figure 25 with a larger value of e 0 which is 0.4. After 1000 relaxation steps, the L 2 norm 
drops 2 orders of magnitude with the larger value of e 0 and a Courant number of 5.0. The 
convergence then stalls at L 2 of approximately 0.0003. A decrease in the Courant number to 
0.5 with the same value of e 0 allows the L 2 norm to drop another 5 orders of magnitude during 
a subsequent 1500 relaxation steps. This behavior is indicative of a low level limit cycle that 
can occur at the larger Courant numbers. 

Figure 27 records the drop in the L 2 norm for the three base region domains for the case 
for a = 0°. The wake region is initialized by inserting the computed values of the solution 
vector from the outflow boundary of the forebody domain into the computational cells in the 
three base region domains. The L 2 norm decreases by more than 4 orders of magnitude in 2600 
relaxation steps. Several Courant numbers were used, as indicated in figure 27. There were no 
limit cycles apparent in these tests, except at a Courant number of 8 x 10®. The eigenvalue 
limiter e 0 during this test was held constant at 0.25. Recall that this value was too small to 
achieve acceptable convergence behavior in the forebody tests. 

The convergence history for the fine grid solution over the forebody at ot = 5 is presented 

in figure 28. The initial condition is interpolated from the initial grid solution. The Courant 
number in this case was greater than or equal to 5.0. The eigenvalue limiter t 0 during this 
test was held constant at 0.2. Convergence in this case stalled at a value of L 2 « 10 5 . This 
level is still nearly 4 orders of magnitude lower than the level that could be obtained in the 
regular grid tests with t a — 0.25 and illustrates the role that truncation error can play in the 
limit cycle behavior. The limit cycle in this case can be even further reduced by reducing the 
Courant number to 0.5, as shown in figure 29, which picks up from the solution state of the 

previous case. ^ y • 

The convergence history for the fine grid solution in the three base domains at a — -5° is 
presented in figure 30. The initial condition is interpolated from the regular grid solution. The 
Courant number in this case was 0.5. The value of the eigenvalue limiter e 0 is 0.01 which is an 
order of magnitude smaller than the value used in the forebody region. Apparently, the greater 
influence of viscous dissipation in the wake and the larger region of supersonic flow surrounding 
the wake core contribute to the improved stability and convergence of the computation in this 
region. 

The two-dimensional inlet calculations did not require any of the parameter adjustments 
used on the blunt body calculations. In particular, there was no “stalling” of convergence 
with large Courant numbers and the required eigenvalue limiter was 2 orders of magnitude 
smaller than used for the blunt body computations. Here again, the absence of any extensive 
stagnation region and the predominance of supersonic flow are believed to improve the stability 
and convergence (fig. 15) of the computation in the inlet relative to the blunt body flow 
computations. 


Concluding Remarks 

An upwind-biased, point- implicit relaxation algorithm for obtaining the numerical solution 
to the governing equations for three-dimensional, viscous, compressible, perfect-gas flows has 
been described. The algorithm is derived by using a finite- volume formulation in which the 
inviscid components of flux across cell walls are described with Roe s averaging and Harten s 
entropy fix with second-order corrections based on Yee s symmetric total variation-diminishing 
scheme. Viscous terms are discretized using central differences. The scheme is second-order 
accurate in computational space except at isolated points where the flux limiter restricts a 
second-order antidissipative correction which could introduce catastrophic instabilities. The 
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scheme is second-order accurate in physical space so long as metric coefficients can be assumed 
to vary continuously from cell to cell. Metric coefficients are defined by the ratio of cell areas 
to cell volumes within the context of a finite-volume scheme. 

The relaxation scheme requires only a single level of storage for the unknowns if time 
accuracy is not required. Time-accurate solutions can be constructed by employing additional 
storage levels to save the results of previous iterations. Because of the point-implicit relaxation 
strategy, the algorithm remains stable at large Courant numbers without the necessity of solving 
large block tridiagonal systems. On machines with large memories, the inverse of the Jacobian 
for the right-hand-side residual vector with respect to the solution vector can be saved over 
large blocks of iterations, causing the algorithm to require very little overhead as compared with 
a standard explicit scheme. These features of the relaxation strategy make the algorithm well 
suited for computers employing either vector or parallel architectures. It is also well suited to 
the numerical solution of the governing equations on unstructured grids, although this option 
has not been developed in the present work. 

Two-dimensional flow problems are used to illustrate convergence properties and sensitivity 
of the solution on numerical parameters. The first test case involves laminar, supersonic flow 
through an inlet. The inflow Mach number is 5.0 and two test Reynolds numbers, based on the 
length of the domain, of 1.14 x 10 7 and 1.14 x 10 6 were used. Convergence of the computed 
solution is evaluated through a grid refinement study. A second test case, involving Mach 14 flow 
over a 15° ramp, is also used to evaluate convergence through a grid refinement study. The 
enhanced numerical dissipation provided by increasing the lower limit on eigenvalue magni- 
tude c is demonstrated. Three different symmetric total variation-diminishing (STVD) flux 
limiters were evaluated in the first test problem. The solution given by the third limiter gives 
sharper pressure jumps and lower skin friction predictions than the first two limiters. The 
predictions using this limiter are, generally, in better agreement with the predictions given by 
three other relaxation schemes. 

The application of separate relaxation factors to the inviscid and viscous contributions 
to the Jacobian within the context of point-implicit relaxation has been found to accelerate 
convergence. The inviscid relaxation factor must be greater than or equal to 1.5 in order for the 
second-order accurate flux-limited schemes to converge. As a converged state is approached, a 
value of 1.5 for the inviscid relaxation factor was found to be the optimum in several numerical 
tests. The optimum value for the viscous relaxation factor was problem dependent in the present 
set of tests, though this behavior may be due to a leading-edge singularity in the boundary 
condition of the inlet problem. The optimum value of the viscous relaxation factor in the blunt 
body problems varied between 1.0 and 0.5, with 0.5 giving the best convergence rate as the 
converged state is approached. 

Several numerical tests were conducted on the Aeroassist Flight Experiment (AFE) con- 
figuration, including the base flow region of a wind-tunnel model. Comparisons of computed 
aerodynamic coefficients and surface pressure with another numerical method and with exper- 
imental data show good agreement. Comparisons of computed heat-transfer coefficients with 
experimental data at three different angles of attack show good agreement on the cone flank, 
but reveal a strong influence of the coordinate singularity at the elliptic nose in the stagnation 
region on the computed heating distributions. Grid refinement in the direction normal to the 
boundary layer and grid restructuring in the vicinity of the axis can be used to improve the 
computed heating distribution in the stagnation region. A four-domain grid was used to dis- 
cretize the space surrounding the forebody and sting of the complete AFE wind-tunnel model. 
Computed contour lines pass smoothly across domain boundaries. Heating distributions on the 
sting are in good agreement with experimental data under the assumption of steady, laminar 
flow. A grid refinement for this case shows only slight truncation error effects on the computed 
sting heating. 
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Appendix A 


Cell Geometry 

Expressions for cell volume and cell wall area in a cell-centered scheme in which the 
dependent variables are at cell centers and the independent variables ( x , y, z ) are at cell corners 
are presented in this appendix. A similar derivation may be found in reference 37. 

Let i,j, k be indices of a cell corner at which point the values of the independent variables 
x,y,z are known. Let I,J,K be indices of a cell center at which point the values of the 
dependent variable q are to be determined. Note, for example, that the computational planes 


represented by I — \ and i are equivalent in this indexing system. 

Let 

dri = rjj+i fc — rjjjt+i = [(ijj — *<;) i + (i/r ; — Vq) j + k ] (Ala) 
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The generic differences, 3f , 3^, 5^, used in equations (Al) are defined as follows: 
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where s is a dummy variable representing the independent variables, x,y,z , and £, 77 , ( are 
computational coordinates in the i,j,k directions, respectively. The cell volume Q and cell wall 
directed area o can now be defined by simple vector relations as follows: 



dr 1 
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vi,j,k = 


y 


(A3c) 

Qi,J,K = 

- dx-j 


+ °I,j,K + °I,J,k) 
3 

(A4) 


where the magnitude of a is the cell wall area, the direction of 3 is normal to the cell wall 
toward increasing I, J, K, respectively, and fl is the cell volume. 

The transformation metrics, such as £ x , £ y , £ 2 , can be expressed as ratios of cell wall areas 
to cell volumes. We would like these expressions to be second-order accurate with respect to 
a face-centered index (i.e., i, J,K) or a cell-centered index ( I,J,K ). Equation (A4) for cell 
volume does not meet this criterion, and use of this expression causes jags in contours running 
across an axis singularity. A better formulation utilizes symmetric averages of differences about 
the cell center (/, J, K) as follows: 

o d?8 ' &,J,K + #i+l,J,K + #IJ,K + fc+l) 

IJ ' K 6 
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Note that the 3 differences are second-order accurate with respect to face centers and that the 
s differences are second-order accurate with respect to cell centers. 

The transformation metrics are defined by 
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In the limit as cell volume goes to zero, using equations (Al) through (A6) gives 
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(A 7a) 
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Equations (A7) are only first-order accurate with respect to any point in the /, J, K cell. This is 
because the a terms are functions of 3 differences which are second-order accurate with respect 
to face centers but first-order accurate with respect to cell centers, and ^i,j,k is a function of 
3 differences which are second-order accurate with respect to cell centers but first-order accurate 
with respect to face centers. Symmetric averages of cell volumes or cell faces can be used to 
achieve second-order accurate metrics at either face centers or cell centers. Several options 
were studied in reference 1. The recommended formulations for face-centered metrics, which 
are required in the evaluation of the viscous dissipation terms across cell faces, are as follows: 
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The unit normal and tangent vectors for the i faces are defined as follows: 
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The unit normal vectors for the j and k cell faces are defined similarly. If the cell wall area is 
equal to zero, extrapolation of n from interior cell walls may be used. 

In the case of the thin-layer Navier-Stokes equations, only the vectors defined by equa- 
tions (A8a), (A8e), and (A8i) are required, depending on the orientation of the coordinate 
system. The dot product of these vectors with the corresponding unit normal to the cell face, 
as used in equations (4.4) through (4.10), can be approximated as follows: 
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where, for example, a iJK = \a iJjK \ is the magnitude (area) of a i JK . 

Another identity that is useful in the programming of the thin-layer Navier-Stokes equations 
involves a geometric relation between the unit normal to a cell face and the gradient of the 
computational coordinate which defines the cell face. This relation is most easily established 
by examining the limiting behavior of equations (A7) and (A9a) as cell volume goes to zero. 
Consider, for example, the gradient of the computational coordinate £ across a computational 
cell face defined by £ = Constant ( i = Constant) given by 
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ZirLE 




Take the dot product of both sides of equation (All) with n to obtain 
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Substitute equation (A12) back into equation (All) and relate the individual vector components 
to show that 



^ n s)i,j,K 
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where s is a dummy variable for x, y, or 2 . Similar relations for the other two coordinate 
directions can be derived in like manner to obtain 
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Appendix B 


Definition of Matrices A, R, A, and R -1 

The Jacobian of g with respect to q is expressed as 
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The similarity transformation matrices R -1 and R, composed of the left and right eigenvectors 
of A, are defined as 
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Note that RR 1 - R J R = I. The diagonal matrix of eigenvalues of A is defined by 
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The matrix definitions presented in equations (Bl) through (B4) need to be evaluated at 
cell interfaces. Consequently, the unit normal and tangential vector components are defined 
with respect to cell walls by using equation (A9). The velocity components in the n, f, and m 
directions are defined as 
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Other parameters required to complete the definition of equations (Bl) through (B4) include 
the sound speed a given by 

a 2 = 2 = (3(H - a) (B6) 

P 

and the kinetic energy per unit mass a given by 

a = ~ (u • u) (B7) 

£ 

The parameter (3 is related to the perfect-gas specific heat ratio 7 by 

/? = 7-l ( B8 ) 

Roe’s averaging is used to define the elements of R, R -1 , and A so that Roe’s Property U 
is maintained. (See ref. 2.) Property U is defined by the relation 

RjA/Rf^qL - qi-i) = El - gi-l ( B9 ) 

where the grid indexing system described in appendix A and averaging in the l direction have 
been used. Equation (B9) is not generally satisfied exactly by standard schemes for obtaining 
average values at the cell interface l based on values at adjacent cell centers L and L - 1, 
because of the nonlinear relation between g and q. Early numerical tests showed that the 
sharpest normal shock captures at hypersonic conditions were obtained using Roe s averaging. 
Roe’s averaged values at cell wall l are defined as follows: 

u ; = c 2 (ciu £ + u^-x) 

Hi = c 2 (ciHl + H l _ 1 ) 

-V 0 

where 

o = ifA 1 ' 2 


(BIO) 
(Bl 1) 

(B12) 

(B13a) 


(B13b) 



Appendix C 


Definition of B 11 for Point-Implicit Treatment of Viscous Terms 

The point-implicit treatment of the governing equations considers only the functional 
dependence of the conservation laws on the cell center which is undergoing a relaxation step. 
Consequently, only derivatives across a cell face contribute to the point-implicit treatment 
because these are the only terms which depend on the cell center. Furthermore, given a 
rectangular ordering of mesh points in which the face normal is directed in the direction of 
increasing index, there is a sign change in the definition of B as compared with B j jr 

because of the different functional dependence of and 011 the cell center as can 

be seen in equations (4.12a) and (4.12b). Following the definitions given in equations (4.8) and 
(4.10) through (4.13), one can define B i i as 


where 


B LL = ~ 


dq l dq L 
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Appendix D 

Definition of Time Step St 

The Courant number is defined as the ratio of a computational time step St to the real, 
physical time required for an arbitrary perturbation in a flow field to traverse a computational 
cell. For a given value of St, the Courant number varies from cell to cell and the largest Courant 
numbers generally occur at the smallest cells. Conversely, if the Courant number is held fixed 
across the domain, then the solution at every cell is advanced at a different computational 
time step. In this case, the temporal evolution of the problem is in error, but the steady-state 
solution is usually independent of the path taken to convergence. 

For most explicit schemes, a Courant number equal to 1.0 represents the largest time 
step which may be used to advance the solution without developing instabilities. Stability 
analyses for explicit schemes are used to define this physically limiting time step. The present 
formulation for time step is modeled after the formulation used in reference 38. Metric 
coefficients in the earlier formulation are replaced by the appropriate ratios of cell volume 
to cell wall area as discussed in appendix A. No stability analysis was performed on the explicit 
counterpart of the present algorithm (as defined by eq. (6.10) with = I). This formulation 
serves only as a point of reference for defining a Courant number. 

The computational time step at cell L ( St)i is given as a function of the Courant number 


C N by 

(St) L c N n L | jRl | + | jR2 | + | i?3 | + fl4 

(Dl) 

where 

^3 

<KS 

r t> 

II 

£ 

(D2a) 


#2 = U L ' V? ft 

(D2b) 


#3 = U L • VCi 

^l = a i(^£z, ■ + • Vr? L + V( L • VCl 

(D2c) 


+ 2|va • V%|+2 |V£l • ^Cl| + 2|V»? l • VC L |) 

(D2d) 

The cell-centered metric coefficients are defined as 



1 + V&+1 ,j,k) 

(D3a) 


VtlL = \ (v m,j,K + Vr)I,j+ i,tf) 

(D3b) 


VCl = ^ (VC/.M + VC/ iM+ i) 

(D3c) 


and the face-centered metrics are defined in equations (A8a), (A8e), and (A8i). The terms R\, 
i? 2 , and i ?3 represent the inverse of the transit time for a convective wave to cross a cell in 
the £, ?/, and £ directions, respectively. The term represents the inverse of the transit time 
for an acoustic wave to cross a characteristic length of the cell. Viscous contributions to the 
formulation of the time step have been ignored. 
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Appendix E 


Boundary Conditions 

Boundary conditions considered in the present work include uniform supersonic inflow, 
supersonic outflow with subsonic boundary layer, no-slip wall, plane of symmetry, and axis 
singularity. Boundary conditions are imposed w ithin the context of a cell-centered scheme 
which utilizes pseudo cells on the opposite side of the boundary. Pseudo cells are not explicitly 
defined, except for the cell face which lies on the real boundary and is shared with the real cell 
inside the computational domain. Pseudo cell volumes are extrapolated from interior cells as 
follows: 


n ° “ nj 


(El) 


Values of dependent variables at pseudo cell centers and values of behind the boundary 
(/ = 0) are required to fully define the problem. When the dependent variables in a pseudo cell 
are a function of the interior cells, the most recent available data from the interior are used. 


Uniform Supersonic Inflow 

Values of dependent variables at the pseudo cell center are fixed at free-stream conditions. 
Values of sq are set equal to 0. 


Supersonic Outflow With Subsonic Boundary Layer 

First-order extrapolation is generally used to define dependent variables at pseudo cell 
centers and is a more stable boundary condition in the highly nonlinear initialization phase of 
the solution process. Second-order extrapolation has been used in some tests, but the results 
were not significantly different from those obtained with the first-onto methods. Positive 
definite quantities are extrapolated to second-order accuracy as in equation (El). First-order 
extrapolation is used to define sq at the pseudo cell wall behind the boundary so that sq = sp 
Because of the upwind nature of the approximation scheme and the hyperbolic/parabolic nature 
of the outflow boundary, the influence of the pseudo cell specification on the interior points is 
very weak. 


No-Slip Wall 

A zero normal pressure gradient is enforced by imposing the pressure at the boundary cell 
onto the reflected pseudo cell. Adiabatic wall boundary conditions are enforced by imposing the 
total enthalpy at the boundary cell onto the reflected pseudo cell. Cold wall boundary conditions 
and no-slip velocity boundary conditions are enforced in one of two ways. The simplest and 
most stable method is to impose the exact wall boundary conditions at the pseudo cell center, 
which is in fact half a cell away from the real boundary. In practice, if the grid is fine enough 
to adequately resolve boundary-layer gradients, then the offset of the pseudo cell center from 
the boundary is small and the influence of this offset is minimal. (Grid resolution is judged 
adequate when the cell height Ah at the wall boundary is approximately equal to the local 
value of / pa.) A numerically more accurate boundary condition is to specify the pseudo cell 
velocity and energy such that the interpolated wall boundary value is correct. Thus, 


u 0 = -Ui (E2) 

e 0 = 2e wall _ e l (E3) 

Equations (E2) and (E3) are more likely to lead to instabilities in the early initialization phase 
of the developing solution. First-order extrapolation is used to define sq at the pseudo cell wall 
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behind the boundary so that so = s*. In effect, this specification of s 0 restricts the flux limiter 
of equation (3.8) to consider only gradients at the surface and in the interior of the domain. 
So long as the min mod function does not return a zero argument, the accuracy of the scheme 
is still formally of second order. 

Symmetry Plane 

Pseudo cells behind the plane of symmetry correspond to actual interior cells. Except 
for the reflection of the velocity component normal to the symmetry plane, the pseudo cells 
are defined by imposing the dependent variables from the corresponding reflected cells in the 
interior. Values of the gradients s 0 are calculated directly from reflected cell center data. 

Axis Singularity 

An axis singularity occurs when the computational domain is constructed by rotating a 
two-dimensional grid around some natural axis of the body or the flow. The flow itself does 
not have to be axisymmetric for this construction to exist. Pie-shaped cells with zero cell-wall 
area on the axis surround the axis singularity. Pseudo cells behind the axis correspond to 
actual interior cells which are reflected across the symmetry plane. Specification of dependent 
variables at pseudo cell centers and gradients at cell walls behind the axis proceed similarly to 
the symmetry plane boundary. However, if J is the index of the boundary cell and JN is the 
number of cells in the semicircle surrounding the axis, then the corresponding pseudo cell is 
a reflection of the interior cell with index JN + 1 — J. Pseudo cell specification has no effect 
on the first-order accurate algorithm (9 = 0) because the common cell wall has zero area. The 
gradient across the axis will influence the outcome of the flux limiter in equation (3.8) for the 
far wall (l = 2) of the boundary cell in the second-order algorithm. 
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Appendix F 

Unstructured Grids 

The formulations contained in this paper are based on a structured grid system (i.e., the 
computational cells are rectangularly ordered with indices i, j, and k, and nearest neighbor 
cells differ in index by one). Unstructured grids are generally composed of triangular cells in 
two dimensions and tetrahedral cells in three dimensions. In the most general case, they 
cannot be ordered as described for structured grids; however, they have an advantage of 
being completely adaptable to very complex geometries and flow structures. In light of recent 
developments in unstructured grid formulations (ref. 39), some comments are offered with 
regard to modifications required in the present analysis for such applications. 

The finite- volume approximation to equation (2.1) for a general, unstructured grid is written 

[^2] +£f m • *™«m = 0 (FI) 

. ot )l 

where <5q = q n+1 - q n and 6t = t n+l - t n . The summation is over all the faces of cell 
L enclosing volume Q. Subscript m indicates a quantity taken on cell face m with surface 
area a m . The quantity n„ t is a unit vector normal to cell face m in a direction facing away 
from the cell center, 

The geometric derivations of metrics in appendix A are based on a rectangularly ordered, 
structured grid system. Some modification of these formulations would be required for general, 
unstructured grids. One point to consider is that the unit normal to a cell face in the 
unstructured grid formulation n m faces away from the cell center, whereas the unit normal 
to a cell face in the structured grid formulation n t faces in the direction of increasing index l. 
This convention accounts for the minus signs appearing in equations (2.4) and (2.5). 

A first-order accurate formulation of the Euler equations on a structured grid is identical 
to the formulation on an unstructured grid. Modifications are required to achieve second-order 
accuracy. Clearly, the present formulation requires an ordered grid system to define sj nn . The 
construction of an equivalent limiter from neighboring cells in an unstructured grid should 
be possible but has not been investigated. In like manner, the formulation of the viscous 
terms, which depends on ordered computational directions in the present case, will need to be 
reformulated in the unstructured grid case. However, the formulations for obtaining second- 
order accuracy with the complete Navier-Stokes equations will only involve modifications to 
the right-hand-side residual vector. The point-implicit relaxation procedure that is defined by 
the formulation of the left- hand-side matrix is independent of grid structure. Consequently, the 
development of point- implicit relaxation will carry over with minimal changes in unstructured 
grid formulations. 
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Appendix G 


Asynchronous Iteration 

A two-dimensional scalar version of the present algorithm was modified to test the feasibility 
of executing on a computer with a massively parallel, asynchronous processing architecture. 
The test is only a simulation in that it was conducted on a serial machine. The environment to 
be simulated is one in which each computational cell has associated with it a single processor 
which is to be kept busy at all times executing its own copy of the master algorithm. When 
the algorithm in a particular cell (processor) requires data from a neighbor cell, it gathers the 
latest available data from local shared memory, regardless of the iteration level of the neighbor 
or where the neighbor is within its own processing of its copy of the algorithm. As the global 
solution develops, individual processors may get out of step. This loss of synchronization 
may arise from variations in instructions at boundary cells, different branches of instructions 
resulting from conditional branches within the algorithm, and potential hardware differences 
in neighboring processors. 

This environment is crudely simulated on a single processor, serial machine by using a 
random number generator to direct the order that cells are relaxed in the computational 
domain. A version of the present algorithm was coded in which a single index L ordered 
the relaxation sweeps through a two-dimensional domain. All aspects of a relaxation step were 
completed at cell L before proceeding to cell L + 1. Convergence histories were recorded for 
the ordered sweeps. The loop over index L was then replaced by a loop over dummy index 
LL, and index L was redefined by a random number generator which ranged from 1 to LMAX, 
where LMAX is the total number of cells. Convergence histories were then recorded for the 
random sweeps, in which case a sweep is defined as one complete pass through the loop. Within 
a given sweep, some cells will likely be called more than one time, and other cells will not be 
called at all. As sweeps continue, the iteration levels of neighbor cells are not generally equal. 
This random offset of iteration levels of neighbor cells executed on a single processor is used to 
simulate the asynchronous iteration of cells on a massively parallel processor. This simulation 
does not model all the possible ways that processors can get out of synch on a real, parallel 
machine. However, it does serve to demonstrate the potential applicability of upwind-biased, 
point-implicit relaxation schemes on asynchronous, massively parallel processors. 

The tests were conducted on two supersonic flows over cylinders. Both tests used initial 
conditions for a body materializing in a uniform, supersonic flow at time t = 0. The first 
test involved Mach 3 inviscid flow and used the first-order scheme with variable time steps 
based on a Courant number of 2.0. The second test involved Mach 1.9 viscous flow (using 
the thin-layer Navier-Stokes approximation) at a Reynolds number of 105 based on cylinder 
radius and used the second-order scheme with variable time steps based on a Courant number 
of 0.75. The convergence histories of the two tests are presented in figure Gl. The random 
sweeps take twice as long as ordered sweeps in the inviscid case and 1.25 times as long in the 
viscous case. In both cases, the converged solutions were identical to at least four significant 
figures, indicating that the order of relaxation contributes to the rate of convergence but not 
to the final converged solution. 

The converged viscous solution and experimental data (ref. 40) for the pressure coefficient 
are presented in figure G2. The comparison in this case is good in the stagnation region but 
breaks down near the outflow boundary. A calculation using the full Navier-Stokes equations 
with slip-flow boundary conditions would be more appropriate for this case, but this capability 
weis not included in the serial version of the code. 

The success of these tests bodes well for the possibility of using an upwind-biased, point- 
implicit relaxation scheme on asynchronous, massively parallel computers. One recognizes that 
there are potential obstacles in the design of a computer that could exploit the synergism 
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(a) Inviscid calculation; Moo = 3.0. (b) Viscous calculation; = 1.9. 

Figure Gl. Convergence history for calculations on serial computer using ordered and random sweeps across 
domain for two-dimensional supersonic flow over cylinder,- 



Figure G2. Experimental data for pressure coefficient and prediction for Moo = 1-9 case with ordered and 
random sweeps. - -- 

of algorithm and architecture described here. Problems involving sharing memory among so 
many processors may overwhelm any potential benefits. However, the real-world challenges 
facing computational fluid dynamicists today, not to mention the computational challenges 
faced by other disciplines, call for algorithms which must deal with millions of computational 
cells that are patched together in ordered subsets or that are completely unstructured. 
Hardware restrictions limit significant speedup of machines with only vector architecture for 
these problems. The optimum blend of vector and parallel features is likely to be algorithm 
dependent, and benchmarking of tomorrow’s supercomputers should not be based strictly on 
today’s algorithms that were optimized for today’s architectures. In anticipation of these 
developments, an upwind-biased, point-implicit relaxation scheme should offer the flexibility 
required for efficient implementation on either vector or massively parallel machines (or some 
combination thereof) utilizing structured or unstructured grids. 
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(b) Skin friction; lower wall. 

Figure 13. Comparison of four algorithms for inlet problem with Re = 1.14 x 10 7 taken from 
reference 32. 
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Figure 15. Convergence history for inlet problem using fine (200 x 200) grid with e 0 = 0.1. 
Bottom of figure shows maximum pressure on upper wall behind reflected shock and 
minimum skin friction (maximum negative quantity) in separated region in front of reflected 
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Figure 18. Aeroassist Flight Experiment (AFE). 
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Figure 19. Experimental data and prediction for pitching moment of aerobrake at Mach 10. 
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(b) a = 5°. 

Figure 20. Prediction and experiment for pressure and heating over AFE at Mach 10. 
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Figure 21. Prediction (laminar) and experimental data for heating on sting supporting AFE 
model in Mach 10 test case. 
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Figure 22. Four-domain grid defining surface and plane of symmetry for AFE test case. 
Domains are approximately divided into forebody, outer wake, shear layer behind shoulder, 
and inner wake core surrounding sting. 
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(a) a - 0°. 



(b) Q = -5°. 

Figure 23. Pressure contours in plane of symmetry about AFE at Mach 10. 


61 
















NASA 

National Aeronautics and 
Space Administration 


Report Documentation Page 


2. Government Accession No. 


1. Report No. 

NASA TP-2953 


4. Title and Subtitle 

An Upwind-Biased, Point-Implicit Relaxation Algorithm for 
Viscous, Compressible Perfect-Gas Flows 


3. Recipient’s Catalog No. 


5. Report Date 

February 1990 

6. Performing Organization Code 


7. Author(s) 

Peter A. Gnoffo 


9. Performing Organization Name and Address 

NASA Langley Research Center 
Hampton, VA 23665-5225 


8. Performing Organization Report No. 

L-16588 

10. Work Unit No. 

506-40-91-02 

11. Contract or Grant No. 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 


13. Type of Report and Period Covered 

Technical Paper 

14. Sponsoring Agency Code 



16. Abstract 

An upwind-biased, point-implicit relaxation algorithm for obtaining the numerical solution to the 
governing equations for three-dimensional, viscous, compressible perfect-gas flows is described. 
The algorithm is derived by using a finite-volume formulation in which the inviscid components 
of flux across cell walls are described with Roe’s averaging and Harten’s entropy fix with second- 
order corrections based on Yee’s symmetric total variation-diminishing scheme. Viscous terms 
are discretized by using central differences. The relaxation strategy is well suited for computers 
employing either vector or parallel architectures. It is also well suited to the numerical solution of 
the governing equations on unstructured grids. Because of the point-implicit relaxation strategy, 
the algorithm remains stable at large Courant numbers without the necessity of solving large, 
block tridiagonal systems. Convergence rates and grid refinement studies are conducted for 
Mach 5 flow through an inlet with a 10° compression ramp and Mach 14 flow over a 15° ramp. 
Predictions for pressure distributions, surface heating, and aerodynamic coefficients compare well 
with experimental data for Mach 10 flow over a blunt body. 


17. Key Words (Suggested by Authors(s)) 

Hypersonic flows 

Aeroassisted space transfer vehicle 
Total variation diminishing 


18. Distribution Statement 

Unclassified — Unlimited 


19. Security Classif. (of this report) 

Unclassified 


20. Security Glassif. (of this page) 

Unclassified 


Subject Category 34 

21. No. of Pages 22. Price 

75 A04 


NASA FORM 1626 OCT 86 

For sale by the National Technical Information Service, Springfield, Virginia 22161*2171 


NASA-Langley, 1990 







